® 


Check for 
updates 


accuracy. The real procedure gauss computes the area under 

the left-hand portion of the normal curve. Algorithm 209 [3] 

may be used for this purpose. If f < 0 or if dfl < lor if df2 <1 

then exit to the label error occurs. 

National Bureau of Standards formulas 26.6.4, 26.6.5, and 
26.6.8 are used for computation of the statistic, and 26.6.15 is 
used for the approximation [2]. 

Thanks to Mary E. Rafter for extensive testing of this proce- 
dure and to the referee for a number of suggestions. 
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begin 

if dfl < 1 V df2 < 1 V f < 0.0 then go to error; 

if f = 0.0 then prob := 1.0 

else 

begin 

real f1, f2, x, ft, vp; 

fl := dfl; f2 := df2; ft := 0.0; 

x := f2/Ff2+flXf); vp := fl + f2 — 2.0; 

if 2 X (dfl+2) = dfl A dfl < mazn then 


begin 
real zz; zz := 1.0 — z; 
for fl := fl — 2.0 step — 2.0 until 1.0 do 
begin 


vp := vp — 2.0; 
ft := xa X vp/fl X (1.0+ft) 


end; 
ft := x Î 0.5xf2) X (1.0+ft) 
end 
else if 2 X (df2 + 2) = df2 A df2 £ mazn then 
begin 
for f2 := f2 — 2.0 step — 2.0 until 1.0 do 
begin 
vp := vp — 2.0; 
fi := x X vp/f2 X (1.0+ft) 
end; 
ft := 1.0 — (1.0—2z) T (0.5XxXf1) X (1.0+f0 
end 
else if dfl + df2 < mazn then 
begin 


real theta, sth, cth, sts, cls, a, b, xt, gamma; 
theta := arctan(sqrt(f1Xf/f2)); 

sth := sin(theia); cth := cos(theta); 

sts := sth Î 2; cts: = cthT 2; 


a := b := 0.0; 
if df2 > 1 then 
begin 


for f2 := f2 — 2.0 step — 2.0 until 2.0 do 
a := cts X (f2—1.0)/f2 X (1.0+a); 
a := sth X cth X (1.0+a) 
end; 
a := theta + a; 
if dfl > 1 then 
begin 
for fi := fl — 2.0 step — 2.0 until 2.0 do 
begin 
vp := vp — 2.0; 
b := sts X vp/fl X (1.0+b) 
end; 
gamma := 1.0; f2 := 0.5 X df2; 


Volume 12 / Number 3 / March, 1969 


for zi := 1.0 step 1.0 until f2 do 
gamma := ri X gamma/(zi—0.5); 
b := gamma X sth X cth Îdf2 X (1.0+b) 
end; 
fi := 1.0 + 0.636619772368 X (b—a); 
comment 0.6366197723675813430755351 -+- = 2.0/r; 
end 
else 
begin 
real cbrf; 
fl := 2.0/(9.0 X fl); f2 := 2.0/(9.0Xf2); 
cbrf := f 1 0.333333333333; 
ft := gauss(— ((1.0—f2)Xcbrf+f1—1.0)/ 
sgrt(f{2Xcbrf T 2+f1)) 
end; 
prob := if ft < 0.0 then 0.0 else ft 
end 
end Fiest 
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procedure SORT(A, i, j); 
value i,j; integer 7, j; 
array A; 

comment This procedure sorts the elements of array A into 
ascending order, so that 


Alk] < Afktl1], k=t,1+1,---,j7-1. 


The method used is similar to QUICKERSORT by R. S. Scowen 
[5], which in turn is similar to an algorithm given by Hibbard 
(2, 3] and to Hoare’s QUICKSORT [4]. QUICKERSORT is used 
as a standard, as it was shown in a recent comparison to be the 
fastest among four ACM algorithms tested [1]. On the Bur- 
roughs B5500 computer, the present algorithm is about 25 
percent faster than QUICKERSORT when tested on ran- 
dom uniform numbers (see Table I) and about 40 percent 


faster on numbers in natural order (1, 2,--- , n), in reverse 
order (n, mn—l,---, 1), and sorted by halves 
(2,4, +--+ n, 1,3, --- ,n—1). QUICKERSORT is slow in sorting 


data with numerous “tied” observations, a problem that can be 
corrected by changing the code to exchange elements afk] > t 
in the lower segment with elements a[g] < ¢ in the upper seg- 
ment. This change gives a better split of the original segment, 
which more than compensates for the additional interchanges. 
In the earlier algorithms, an element with value ¢ was selected 
from the array. Then the array was split into a lower segment 
with all values less than or equal to t and an upper segment with 
all values greater than or equal to ¢, separated by a third seg- 
ment of length one and value t. The method was then applied 
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TABLE I. Sorting Times IN SEcoNDS For SORT AND 
QUICKERSORT, on THE Burrovexus B5500 
CoMPUTER—AVERAGE OF Five TRIALS 
Algorithm 
QUICKERSORT 


Original order and number of items SORT 


Random uniform: 


500 0.48 0.63 
1000 1.02 1.40 
Natural order: 
500 0.29 0.48 
1000 0.62 1.00 
Reverse order: 
500 0.30 0.51 
1000 0.63 1.08 
Sorted by halves: 
500 0.73 1.15 
1000 1.72 2.89 
Constant value: 
500 0.43 10.60 
1000 0.97 41.65 


recursively to the lower and upper segments, continuing until 
all segments were of length one and the data were sorted. The 
present method differs slightly—the middle segment is usually 
missing—since the comparison element with value ¢ is not re- 
moved from the array while splitting. A more important differ- 
ence is that the median of the values of A(z], A[(¢+7)+2], and 
A[j] is used for t, yielding a better estimate of the median value 
for the segment than the single element used in the earlier 
algorithms. Then while searching for a pair of elements to 
exchange, the previously sorted data (initially, A[z]<t< Al[j]) 
are used to bound the search, and the index values are compared 
only when an exchange is about to be made. This leads to a small 
amount of overshoot in the search, adding to the fixed cost of 
splitting a segment but lowering the variable cost. The longest 
segment remaining after splitting a segment o` n has length 
less than or equal to n — 2, rather than n — 1 as in 
QUICKERSORT. 

For efficiency, the upper and lower segments after splitting 
should be of nearly equal length. Thus ¢ should be close to the 
median of the data in the segment to be split. For good statis- 
tical properties, the median estimate should be based on an odd 
number of observations. Three gives an improvement over one 
and the extra effort involved in using five or more observations 
may be worthwhile on long segments, particularly in the early 
stages of a sort. 

Hibbard [3] suggests using an alternative method, such as 
Shell’s [6], to complete the sort on short sequences. An experi- 
mental investigation of this idea using the splitting algorithm 
adopted here showed no improvement in going beyond the final 
stage of Shell’s algorithm, i.e. the familiar “sinking’’ method of 
sorting by interchange of adjacent pairs. The minimum time 
was obtained by sorting sequences of 11 or fewer items by this 
method. Again the number of comparisons is reduced by using 
the data themselves to bound the downward search. This 
requires 


Ali-1] < Ak), tk <j. 


Thus the initial segment cannot be sorted in this way. The 
initial segment is treated as a special case and sorted by the 
splitting algorithm. Because of this feature, the present al- 
gorithm lacks the pure recursive structure of the earlier al- 
gorithms. A 

For n elements to be sorted, where 2% < n < X+, a maximum 
of k elements each are needed in arrays IL and IU. On the B5500 
computer, single-dimensional arrays have a maximum length 
of 1023. Thus the array bounds [0:8] suffice. 
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This algorithm was developed as a ForTRAN subroutine, then 
translated to ALGoL. The original ForTRAN version follows: 


SUBROUTINE SORTIA, II JJ) 
SORTS ARRAY å INTN INCREASING ORDER, FROM ACIT) TO ACJIJ) 
ORDERING IS BY INTEGER SUBTRACTION, THUS FLOATING POINT 

NUMBERS MUST BE IN NORMALIZED FORM. 

ARRAYS TUCK) AND ILIK) PERMIT SORTING UP TO 2+#(K+1)-1 ELEMENTS 
DIMENSION A(1),TUL1L6)eIL (16) 
INTEGER AyT,TT 
M=1 
Isl] 

J=J4 

5 IF(I GE. J) GG TO 70 

10 K=I 
IJ=(J+1}/2 
T=ACTJ) 

IF{ACT) «LE. T) GO TO 20 

ACIJI=ACT) 

A(T)=T 

T=A(TS} 

20 L=J 
IFCACJ) «GE. T} GO TO 40 
ALTJ) =A J) 

A(J)=T 

T=A(1J) 

IF(ACT) «LE. T) GO TO 40 

ACT J=A0T) 

AUE)=T 

T=ACT JI) 

GO TO 40 
20 AC(LI=ACK) 

A(K)=TT 
40 L=L-1 

IF(ACL) 2GT. T) GN TN 40 

TT=A(L) 

SO K=K+1 
IF(AC(K) «LT. T) GO TO 50 
IF(K LE. L) GO TO 30 
IF({L-] .LE. J-K) GO TO 60 
IL(M)=1 
TUIM) =L 
=k 
M=M41 
GG TO 80 

60 IL(M)=K 
TUCM)=J 
Jet 
M=M+1 
GO TN 80 

70 M=M-1 
TF(M .EQ. 0) RETURN 
T=TLOM) 

J=TUCM) 

80 IF(J-I «GE. 11) GO TO 19 
IF({I EQ. II) GOTO 5 
T=I-1 

90 I=T+1 
IFI EQ. J) GO TO 70 
T=A(T+1) 

IF(ACI) LE. T) GO TO 90 

K=] 

100 A(K+1)=A(K) 

K=K-1 

IFT eLTe A(K)) GO TO 100 

AIK+1)=T 

Go TO 90 

END 


This Forrran subroutine was tested on a CDC 6400 computer. 
For random uniform numbers, sorting times divided by n log: n 
were nearly constant at 20.2 X 10-* for 100 < n < 10,000, with 
a time of 0.202 seconds for 1000 items. This subroutine was also 
hand-compiled for the same computer to produce a more efficient 
machine code. In this version the constant of proportionality 
was 5.2 X 10-6, with a time of 0.052 seconds for 1000 items. In 
both cases, integer comparisons were used to order normalized 
floating-point numbers. 
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begin 
real {, tt; 
integer ti, ij, k, L, m; 
integer array IL, IU[0:8]; 
m := 0; ii := i; go to l4; 
Li: ij := (+j) + 2; t:= Aftj]; k:=i; L:=j; 
if Ali] > i then 
begin A[ij] := 
if A[j] < t then 
begin 
Alij] := Alj]; Aljl:= t; t:= 
if A[i] > t then 
begin A{ij] := 
end; 
L2: L:= L -— 1; 
if A[L} > t then go to L2; 
= A[Z]; 
L: k:=k+1; 
if A[k] < t then go to L3; 
if k < L then 
begin A[L] := Afk]; Al[k] := tt; go to L2 end; 
if L — i > j — k then 
begin ILim] := i; IU[m] := L; i := kend 
else 
begin ILim] := k; IU[m] := j; j := Lend; 
m:=m+1; i 
I4: ifj — i > 10 then go to Ll; 
if i = ii then 
begin if i < j then go to L1 end; 
for 7 := 7+ 1 step 1 until j do 


Alil; Ali} := = A[zj] end; 


Alig]; 


Aft]; Alt] := t; t:= Al[tj] end 


begin 
t:= Ali]; k:=i-1; 
if A[k] > t then 
begin 
L5: A[k+1] := Alk]; k:i=k—1; 

if A[k] > t then go to L5; 
A[k+1] := t 
end 

end; 

m := m — 1; ifm > 0 then 

- begin ¿i := IL[m]; j := IU[m]; go to L4 end- 


end SORT 


REMARK ON ALGORITHM 329 [G6] 
DISTRIBUTION OF INDISTINGUISHABLE. OB- 
JECTS INTO DISTINGUISHABLE SLOTS [Robert 
R. Fenichel, Comm. ACM 11 (June 1968), 430] 
M. Gray (Recd. 20 Sept. 1968) 
Computing Science Department, University of Adelaide, 
South Australia 
As the procedure stands it is incorrect. Preceding 
end skip 99,189,198, ete. 
the following statement should be inserted: 
if qik] ~ 0 then LefimostZero := k +1 
Thus the compound statement becomes: 
begin 
LeftmostZero := LeftmostZero —1; 
qlk] := g[LeftmostZero] — 1; 
q(LefimostZero] := 0; 
qlLeftmostZero—1] := g{LeftmostZero—1] + 1; 
if gik] = 0 then LefimosiZero := k +1 
end skip 99, 189, 198, etc. 
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REMARK ON ALGORITH 339 [C6] 

AN ALGOL PROCEDURE FOR THE FAST FOURIER 

TRANSFORM WITH ARBITRARY FACTORS 
{Richard C. Singleton, Comm. ACM 11 (Nov. 1968), 
776] 

Ricwarp C. SINGLETON (Recd. 27 Nov. 1968) 

Stanford Research Institute, Menlo Park, CA 94025 


KEY WORDS AND PHRASES: fast Fourier transform, complex 
Fourier transform, multivariate Fourier transform, Fourier 
series, harmonic analysis, spectral analysis, orthogonal poly- 
nomials, orthogonal transformation, virtual core memory, 
permutation 
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On page 778, column 2, the 7th and 6th lines from the bottom 
should be corrected to read: 

LJ: jj := Cli—2] + jj; if jj = C[i—1] then 

begin? := i — 1; jj := jj — C[i]; go to LJ end; 
On page 779, column 1, the 9th and 8th lines from the bottom 
should be corrected to read: 

LX: jj := Dik+1) + jj; if jj 2 Dik] then 

begin jj := jj — Dik}; k:= k +1; goto LX end; 
In both cases jj was originally used as the controlled variable of 
a for clause and thus was undefined after exit; the corrections 
preserve the value of jj for later use. 

If the user prefers to compute constants with library functions, 
line 5 in column 2 on page 777 may be replaced by: 

rad := 8.0 X arctan(1.0); c30 := sgrt(0.75); 

Algorithms 338 (Comm. ACM 11 (Nov. 1968), 773] and 339 were 
punched from the printed page and tested on the CDC 6400 
ALGOL compiler. After changing a colon to a semicolon at the end 
of line 37 in column 2 on page 775, the test results agreed with 
those obtained earlier with this compiler. 

When computing a single-variate Fourier transform of real 
data, procedure REALTRAN may be used with procedure FFT 
(Algorithm 339) to reduce computing time. Two versions of 
REALTRAN have been given (Algorithms 338 and 345 [Comm. 
ACM 12 (Mar. 1969), 179-184]); the first version is the faster of 
the two, but the second should be used if arithmetic results for 
real quantities are truncated rather than rounded. 

In describing the evaluation of a real Fourier series, in the 
middle of column 2 on page 776, the necessary steps of reversing 
the signs of the B array values both before and after calling FFT 
were omitted. The correct steps, including scaling, are as follows: 

REALTRAN(A B, n, true); 

for j := n — 1 step —1 until 0 do Bj] := 

FFT (A, B, n, n, n); 

forj := n — 1 step —1 until 0 do 

begin A[j] := 0.5 X A[j]; BU] := 


—B[j]; 


—0.5 X Bij] end; 


The policy concerning the contributions of algorithms to 
Communications of the ACM appears, most recently, in the 
January 1969 issue, page 39. A contribution should be in the 
form of an algorithm, a certification, or a remark. An al- 
gorithm must normally be written in the ALGOL 60 Refer- 
ence Language or in USASI Standard FORTRAN or Basic 
FORTRAN. 
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